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Abstract 

We discuss a numerical method for convection-difTusion-reaction problems with a free 
boundary in ID. The method is based on the numerical modelling of the interface evolution, 
the transformation to a fixed domain problem and the approximation by an ODE system. The 
interface evolution is modelled by means of the local shape of the corresponding travelling 
wave solution. The method can be applied to many free boundary problems with a finite 
speed of the interface. The presented method can also approximate some problems with an 
infinite speed of the interface for damped travelling wave type solutions. In the numerical 
experiments we compare our numerical solution with the analytical ones for some problems. 

Keywords: interface modelling, free boundary problems, porous media equations, degenerate 
diffusion-adsorption 

1 Introduction 

We shall consider convection-diffusion-reaction equations of the type 

dtu = dlu" + bod^u' + Cqu™ (1) 

on the halfline x e (0, oo) for t > 0, with the boundary conditions 

a) u = 0(<), or b) - d^u" + bov? = q{t). (2) 

at the boundary a; = 0. Here, 5o and cq are constants. The initial condition is of the form 

u{x,Q)^Ufi{x)>Q, (3) 

where uq{x) > for x G (0, Lq) and uq(x) — for x > Lq. We are only looking for nonnegative 
solutions since negative solutions have no physical interpretation. We assume that the parameters 
n, 7 and m are nonnegative real numbers and guarantee the existence of a function s{t) such that 
the solution is positive in (0,s(i)) and u{x,t) — for x > s{t). The function s{t) represents the 
time evolution of the interface. Such solutions have been intensively studied in the last 30-years 
in many papers, e.g., O.A. Olejnik, A.S. Kalashnikov, Czou-Ju-Lin [22, G.J.Barenblatt [3], D.G. 
Aronson |T], B.H.Gilding [TI1[T2], L.A.Peletier [55^, R.Kersner [H], B.H.Gilding and R.Kersner 
[12 E] etc. (see also the list of references in the previous papers). 

We shall discuss the numerical approximation for ([l])-([3]) based on interface modelling. In these 
problems the free boundary s(t) is not given explicitly in the model and its determination is 
included in the solution. Solutions to the linear convection-diffusion-adsorption problems (n — 
TO = 7 = 1) in ([l]) do not show the interface property. There, an initially localized support of 
the initial profile spreads with an infinite speed (s(t) = oo). Degeneration of the elliptic part 
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(n > 1) is the main reason for the appearance of the interface. However, also in the regular case 
(n = 1) the interface will appear under special properties of the adsorption represented by cq and 
TO. Very substantial results (necessary and sufficient conditions for n, m, 7 in Q) to guarantee the 
appearance of the interface have been obtained in many papers, e.g. B.Song ([5S]), B.H. Gilding and 
R.Kersner ([T51 E]), see also list of references there. If 60 7^ 0, cq < and m, n, 7 > the existence 
of an interface is guaranteed when min(n,7) > min(TO, 1). This condition is also necessary (see 
[2S]). If 5o = and Cq = in ([T|) then the condition 71 > 1 is the real criterion whether or not 
the interface appears. An equation with this property is called "porous media equation" and 
an analytical solution has been found by Barenblatt [3- . This equation is a mathematical model 
for adiabatic infiltration of a gas into a porous media. Its solution is not a classical one (but a 
generalized, weak solution, defined by a corresponding integral identity) since dxU is not continuous 
along the interface, i.e. x = s{t). The Barenblatt-Pattle solution of the Cauchy problem for the 
porous media equation 

dtu = dlu'\ n > 1, x e (-cx),+oo) (4) 



is of the form 



where 



2n(n+l) 
n — 1 

and the corresponding initial condition is 

0, for|x|>s(0). 

This analytical solution stimulated the theory of existence and uniqueness for the more general 
"porous media" type problems. More detailed information on this topic can be found in |13l 114] 
and the references therein. 

The solution of ^ for n = 6 and t e (0, 200) is drawn in Fig.l (in 10 equidistant time sections) and 
the evolution of its interface is drawn in Fig. 2. We can readily verify the validity of the equation 

m = — ^a,u"-iu^,(,), (6) 

n — 1 ^ w 

which is the governing ODE model for this iterface. Solutions of the problem (l)-(3) have the 
shape of damped travelling waves with eventually sharp fronts at the interface. This makes the 
numerical approximation very difficult (the movement of the interface is not priori known) . There 
exist numerical methods for free boundary problems that estimat the position of the interface and 
then apply a local adaptation of the space discretization, which are very expensive but have a 
limited precision (see e.g. |T71 [TS] etc.). Some other results related to our method can be found 
in HE]. 



The main goal of this paper is the construction of an efficient and precise numerical method 
for solving "porous media" type problems ([l])-(|3]) . This method is based on the modelling (in 
terms of an ODE) of the time evolution of the interface s{t) , which significantly depends on the 
local profile of the solution near the interface. The unknown s will be explicitly included into the 
solution of the problem. Then, we shall look for the solution in the interval (0, s{t)) only (instead 
of (0,i), where L > Lq). Consequently, we will transform ([T]) to the fixed domain (0, 1) (using 
y — ) and invoke a special space discretization (finer in the neighbourhood of the point y = 1 
which corresponds to the interface x = s{t)). This approximation of ([l]) results in a corresponding 
system of ODEs including the uknown s{t). This system is completed with the ODE model for the 
interface (of the same type as ^ for Q). The final ODE system is solved with an appropriate 
solver for stiff ODE systems. Thus, the modelling of the interface evolution is the crucial point 
for the construction of our numerical method. In this way we can approximate the sharp front of 
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the solution at the interface with higher accuracy. In our numerical experiments we have used the 
solver 'odel5s' from MATLAB hbrary. 

In Section 1 we describe the method in more detail. In Section 2 we construct the interface model 
for ([T|) under some assumptions on parameters n, m and 7. In section 4 we construct the numerical 
approximation for ([ij . In Section 5 we extend the results from the previous sections to the more 
general PDE 

dtu = dla{u) + d^b{u) + c{u) (7) 

satisfying a(0) = 6(0) = c(0) — 0, a'(0) = and preserving similar interface property of the 
solution as ([l]) . Also more general mathematical models involving a{t,u),b{t,x,u),c{t,x,u) can 
be included too. These types of equations describe a large part of convection-difFusion-adsorption 
problems with many important practical applications. In Section 6 we discuss the interface mod- 
elling in the non degenerate case a'(0) > 0. These equations are of the "oxygen-consumption" 
type, where appearence of the interface is generated by strong adsorption represented by c(u). In 
Section 7 we present some numerical experiments and we present comparisons with the analytical 
solutions. 



2 Idea of the numerical approximation 

For a good numerical approximation of ([ij it would be desirable to use space discretization only in 
the moving domain (0, s(t)) and to use also moving grid points which follow the interface and which 
are more dense at the sharp fronts near the interface - see Fig.l. For this purpose we formally 
consider (jll in the domain (0,.s(t)) which by Landau's transformation, y — is transformed it 
to the fixed domain (0, 1). Then, for u{y,t) := u{x,t) we obtain the equation 

_ g 

ut = C{u) - -ydyu, y e (0, 1) (8) 

where C{u) is the rhs of M and C{u) is its transform in the y-variable. If we can model s in 
terms of dyu{l,t), resp. d^{l,t) (which describe local properties of u in a 'small' neighbourhood 
of s{t)), e.g., in the form (see (jlj) 

s = B{dyu{l,t)), resp. s = B{dyu{l,t),d^u{l,t)), (9) 

then we can eliminate s in ([s]) by means of the rhs of (|9| . Then the complete system consists from 
([s]) and We approximate it using the space discretization {yi}fLi S (0, 1) of the y- variable. 
The resulting system of ODEs is stiff and we solve it by a corresponding ODE solver. Now, the 
unknown interface is explicitly included in the solution. The grid points {j/ijfli correspond to the 
moving grid points {xi{t)}fL^ g {0,s{t)), with x^it) = s{t) for ([l]). 

The method used cannot be extended to more space dimensions (except of radial symmetric 
ones) . The mathematical model for the interface movement (its normal velocity) depends not only 
on local properties of the solution (e.g., the normal space derivative) but also on the geometry of 
the support of the solution, which in ID is very simple. 

3 Modelling of the interface 

The appearance of the interface in ([T|) is closely related to the existence of the travelling wave 
solution (see [T31II1] ) of the form u{x, t) = f{^), — x — at where 

/(e)=0, fore>6 
/(e)>0, fore<^o 

for some ^0 G (~oo, -l-oo). The function / has to be determined from the ODE 

-af = n/"- V" + nin - l)r-\f? + bnlP-'if))' + Cof, C e (-00, 00). 
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The function 9 defined as 

oum = -irno, 

has to satisfy the nonhnear Volterra equation (see ^3J ) 



eis) ^as + boS-< - nco / ^^^^ dr. (10) 



In [T5] (Lemma 13) one proves: There exists a solution of ( 10 ) for all a > a* (for suitable a* > 0) 
provided the conditions (i)-(iii) below with parameters n > 1, m > —n and 7 > 1 are fulfilled. 
Moreover, there exists a a** > a* so that 

d{s) w Bqs'i as s \ (11) 

for some and g > depending on n^m,^. These conditions are 

(i) Co < 0; (in this case (jf = min{(n + m)/2, 1}) 

{a) Co — 0; (in this case 9=1) 

{Hi) Co > 0, m + n > 2; (in this case q — 1). 



The conditions (i)-(iii) are also necessary for the existence of the solution of (10). 

The case < n < 1 is also discussed in [13j (Lemma 13), but we cannot aply it in our analysis. 

The relation ( [TT] ) is equivalent to 

e{s) = 6los« + o(s9), with s-«o(s«) -> 0, for s \ 0. 



As a consequence of (111, also the front of the semi- wave / (localized at the point ^0), satisfies 
the approximation (see Remark 1 below) 

/(O « d^(?o - 6^ fore /Co. (12) 

We shall use this property in the modelling of s{t) in terms of local properties of the solution 
at the interface (e.g. dxu(s(t),t)) and model parameters. This will be substantially used in our 
numerical method. In fact, we shall assume the asymptotic property (12 1 to hold also for the first 
and second derivatives of / (see ([T3| below). 



Remark 1. From (jIT|and 0{f{S^)) := ^/"(e)) (provided / is smooth, decreasing and = 
for £, > ^0 - "travelling semi- wave") it follows 

n 



and consequently, by integration over (^,^0)1 we obtain (see ( |12[ |), 

/(C)-[^5^^m-0^ where 

n n — q n 

under the assumption n ^ q. 

The development of an interface model of the type (|9]) in a rigorous mathematical way requires 
a very deep analysis. Moreover, there are no results available in this respect (except of a few, very 
special cases). We will construct such a model under very strong smoothness assumptions on the 
solution in the interval (0, s{t)] (up to the free boundary) and the assumption ( 12 1 concerning the 
shape of travelling semi wave corresponding to (1). In fact we assume (12 1 in its stronger form 

/(O « d^i^o - 0^ no « 0'{^o - 0^-' for e / ^0, 

[r]'(0 « d"'3n/3(eo - 6"^"' fore/eo, (13) 
[/"]"(0 « d-^ndinP - mo - 0''^-^ for C / Co 
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where f3 and d are unknown parameters to be determined. 



Our main result in this section is 



Theorem 1. Suppose (13 1 and n> 1, 0<m<l,7>l. When the solution of (1) is regular up 



to the interface (the equation is satisfied in the limit sense at the interface), then the interface s 
of (1) is governed by ODE 



s{t) = Q{w) 



(n-l) 



dxW 
dxW 
dxW 



Co 



(n-l) 
n 

'(n-l) 



G(7)6o 
G(7)&o 
G(7)&o 



for m + 
for m 4 
if Co 



n = 2; 
n>2 
= 



sit) 



(14) 



where j3 — -^^^ u — and G{r) = for r > 1, G(l) 
Proof. Our starting point is the condition u{s{t),t) 
equation with respect to t we obtain formally 



for all t e (0,r). Differentiating this 



dtu . 



dx 



in the sense x /* s(t), 



(15) 



which cannot be used practically, since usually we have dxU — > — oo for x Z' s(i). But invoking 
the assumption (13 1, we can conclude that dxU^ is finite for x Z' s(i) for some /3 > . Using the 
smoot hnes s of the solution (the equation (1) holds up to the interface x — s{t)) we can eliminate 
dtu in ( 15 1 by the rhs of (1) and then we insert u — d^{s{t) — x)^ + oq in the small neighbourhood 
of the interface (for x /' s{t)). Due to (13 1, we obtain 



s(t) = lim 



1 



ys #(s-2;)/3-l + 



Ol 



-{d"'^n/3(n/3- l)(s(<) - a;) 



n(3~2 



(16) 



bnd^f'-f/Sis ~ x)^'3-i ^ cod™'^(s - x)""'' + 02 + 03 + 04} 



with the lower order terms represented by 01—04 ( because of (13l), where 



Oo 



0{{S-Xf), 0^=0iis~xf-'), 02 = 0((s-x)"'3-2), 



03 = o((5-a;)''^-i), 04 = o((s - a;)"'^) 



m/3 \ 



The assumtions on n,m and 7 guarantee the existence of the interface, i.e. |s(t)| < 00 because 
of [13] (Lemma 13). Our goal is to find such a /? > that this will be achieved. There are 
more combinations of parameters n, m, 7 to guarantee that such a /? could be found. Deviding 
the nominator on the rhs in (16 1 by {s — x)^~^ , we obtain three terms containing s — x with 
the exponents (n — l)/3 — 1; (7 — l)/3; (m — 1)(3 + 1 (corresponing to diffusion, convection and 
reactive part) . To guarantee \s{t)\ < 00, we have to choose f3 > such that all three exponents 
are nonnegative and at least one of them has to be zero. The case /3 < doesn't correspond to 
the required property of the semiwave. Then, for n + m > 2 the choice (3 = fulfiUes our 
requirements. This choice also includes the case cq = . Then from ( 16 1 we deduce the formula 
for s{t) which contains the uknown parameter d -see ( |12[ ). In the following we shall determine d . 
For this purpose consider the new unknown w defined by u = and rewrite (1) in terms of w. 
In the small neighbourhood of the interface {x s{t)) we have w{x,t) « d{s{t) — x) due to (13 1 
and hence 

. w{x, t) — w{s{t), t) 



s{t) 



d\^0, as X Z s(t). 



Prom this we conclude that dxW —d when x /" s{t). 

Now, we substitute d — —dxW into (16 1 and we obtain the required ODE model for the interface. 



Remark 2. It is very difficult to verify the assumptions (13 1. However consequences from 



them, leading to the interface models, are verified in some examples in Section 7, where the an- 



alytical solution is available. Since assumptions (13 1 could be violated in some model problems 
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our modelling of the interface should be justified numerically in practical applications. 
Remark 3. The analysis used in choosing (3 so that \s{t)\ < oo is rough and cannot distinguish 
the sign of 60 and cq which is important for the existence of the interface. The analysis for the 
adsorption and reaction is the same. From (i) and (iii) it follows that the interface exists also in 
the case n + m < 2 provided cq < (i.e. in the adsorption case), while for cq > the interface 
doesn't exist. 

Remark 4. In the case n + m > 2 it follows from (i),(iii) that q = 1. Then, our /3 obtained 
in the proof of Theorem 1 coincides with the one from Remark 1 what supports our arguments. 
In fact, the assumptions in Theorem 1 include also the additional solution /3 = jtz^ in the case 
n + m > 2. However, this /3 is excluded by the fact q — I (see (i), (iii) and Remark 1). It means 
that our analysis leading to Theorem 1 is not satisfactory for the unique determination of f3. 



4 Numerical approximation of (1) 

The mathematical model for the time evolution of s given by Theorem 1 is of practical use, 
since dxW is finite and nonzero. From this reason we also rewrite (1) in terms of w using the 
transformation u — where B = . Then, we obtain 



(17) 



We have to solve this equation in the moving domain (0,s(t)). We transform (17 1 to the fixed 
domain (0,1) using Landau's transformation y — with w(y,t) — w{x,t). For the sake of 
simplicity, in the following we write w in the place of w. Then, we obtain 



dtW = n{n(3 — 1) 



w 



(n-l)/3 



-dyW 



w 



(n-l)/3-l 



,(7-l)/3 



-dyW- 



This equation has to be completed with the ODE for s given by Theorem 1 and s at the rhs has 



to be eliminated by means of the rhs of ( 14 1 . Then, we have to solve the coupled PDE and ODE 
system 

^(n-l)/3-l y;(7-l)/3 

dtW ^ n{n(3 - I) ^ — dyW + n (dyw) +6o7- 



-dyW~ 



s = Q{w) 



(18) 
(19) 



on a fixed domain y £ (0, 1), t > 0, where Q is from (14 1 (there, d^w has to be replaced by '^^) 



The PDE (18 1 at the point {y,t) depends also on the point (l,t) because of Q{w). Consequently, 
it is non-local. In the numerical approximation of (18l-([l9| we use space discretization (method 
of lines) which results in an ODE system. We generally use a nonequidistant partition of (0, 1). 
Let ai > 0, Vi = 1,...,-/V, ao ^ and consider the grid points yi — X]}=o'-'^j' ~ 0,...,N. 
We approximate ■w{yi,t) Ci{t). Denote by £i{y) the Lagrange polynomial of the second order 
through the points (yi„i, Ci_i), (y^, C^) and (j/i+i, C^+i). We approximate dyw{yi, t) and dyw{yi, t) 

by the corresponding expressions given by ^^\y=yi and ^\y=yi- 
We approximate the derivative dyw{l,t) appearing in Q by 

Mb, 



dyW{l,t) 



dy 
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where Ib is the second order Lagrange polynomial throudh the points {yN-2, Cn~2) and {i/n-i, Ctv-i), 
(1, 0). Our ODE system now reads 



Co Q{C) dii 

s — Q{C), where dyw{l,t) ^ ^^^\y^i (21) 

ay 

for i = 1, ...,iV-l, with /? = In the case of Dirichlct boundary conditions this is the resulting 
system of ODE, where we use Co(t) = ^(i) (see (|2|, case a)). In the case of a Robin condition (see 
([2]), case b)), the system (p0|-(21 1 has to by completed by the additional ODE at the point y = 0- 



In this case we make use of a ghost point y^i = — ai. The missing value C_i is determined from 
the boundary condition (|2|, (case b), by means of q and Ci from the relation 







Then, the additional ODE in y = is the same as (20 1 for i = using C_i. 

In some physical problems it is important to know the support of the solution which in ID is 
given by the position of the interface. This is explicitly included in our approximation represented 



by (20 1- (21 1. This gives us the possibility to determine s{t) with higher accuracy then the 



approximations without interface modelling. 



Remark 5. The system (20l-(21 1 is stiff and a corresponding ODE solver has to be used. This 
system is a good starting point to solve the original free boundary value problem. The first 
advantage is that we approximate the solution on its support only. The second one lies in the 
space discretization. In fact, we have moving grid points (in the original x- variable) by means of 
which sharp fronts near the interface can be approximated very well for all t. We shall demonstrate 
this in our numerical experiments presented in Section 7 . The presented numerical approximation 
requires positive initial profile on its support (0, s(0)) . If our initial profile is zero (provided nonzero 
boundary condition at a; = is considered), then we have to start with the "small" positive 
auxiliary initial profile on the " small " domain (0,s(0)). Also some compatibility between the 
initial and boundary conditions (at the point (t, x) = (0, s(0))) help to avoid numerical instabilities 
in the ODE solver at the starting time point. This "small" regularization does not change the 
expected solution significantly. The efficiency of the numerical method considerably depends on 



the strategy of grid points distribution. This allows reducing the size of the system (20 1 without 
decreasing the accuracy. 

The used space discretization of the PDF and its approximation it by an ODE- system is of the 
second order and the order of ODE solver is regulated by its tolerance parameters. 



5 Extension of the method to PDE ([T]) 

From the theory of porous media type solutions one can notice that the evolution of the front of the 
solution substantially depends on the order of degeneration and adsorption and on local properties 
of the solution. Therefore, the same idea of the numerical approximation used on ([T|) can also be 
applied in a more general case of ([t]), where we replace a{s) a{t,s), b{s) b{x,t,s), c(s) <-^- 
c{x,t,s), provided a,b and c are asymptotically polynomial with respect to s at the point s — 0. 
More precisely, we assume 

a{t,u) = g{t,u)u^ , b{x,t,u) — bo{x,t)u'^ , c{x,t,u) — CQ{x,t)p{u)u^ (22) 

with git, 0) 7^ 0, p(0) 7^ and g G 



7 



Modelling the interface in this setting (generahzation of (|7|), we proceed similarly as in Section 2 
and obtain. 

Theorem 2. Let the assumptions of Theorem 1 hold true. Moreover, suppose (22 1. Then, the 
governing ODE for the interface in the genaralized case of ^ is given by 

f -^a,u;+^5^i^f^-G(7)feo(5W,i) ,form + n = 2; 
s(t) = OH = < -j^d,w-Gi^)boisit),t) ,form + n>2, (23) 

[ -j^d,w-G{j)bo{s{t),t) ,forco = 0, 

where u — with (3 = and G{r) = for r > 1, G(l) = 1 and u = with f3 = 

Using this ODE model for the interface, we approximate the considered PDE similarly as in the 

case of ([T]). 



6 Modelling of the interface in the nondegenerated case 

dua{t,0) > 0. 

Cosider the generalization of ([t]) as in the previous section with the nondegenerate elliptic part 
dua{t, 0) > 0. As it was mentioned in the introduction, the interface can appear also in the 
nondegenerated case when adsorption is sufficiently strong for small values of u . As an example 
we can take oxygen-consumption type problem 

dtU = Ddlu-H{u), x£{0,oo),t>0 (24) 

with u{x, 0) = uo{x) > for x e (0, Lq), uo(x) = for > Lq, where D is the diffusion coefficient 
and 

1 for 2 > 



H{z) = 



for 2 



represents the consumption of the oxygen with concentration u. This model problem has also an 
interface with finite speed and has been analysed some decades ago (see, e.g., [S]). There exists 
an analytical solution wich is a polynomial of the second order at the front (near the interface). 
This motivates us to loolc for a polynomial solution at the interface also in the more general case 
([t]) considered in the previous section. The order of this polynomial can be determined from the 
information that the considered problem has a finite speed of interface propagation. Then, for the 
interface modelling we proceed similarly as in the degenerated case. The existence of the interface 
(in case of ([t])) has been established in [T^ (Theorem 12, Assertion (i)). 

For simplicity we consider a generalization of the mathematical model ([t]) (of oxygen-diffusion 
type) under the assumptions 

a~ao{x,t)u, b — bo{x,t)u, c = CQ{x,t)H{u)u"^ (25) 

with aQ{x,t) > (5 > 0, ca{x,t) < -5, < m < 1. 

A more special case when c — cqH{u) has been treated numerically in [7] and we have extend this 
idea to the more general case of ^ , where ( 25 ) is fulfilled. Another approximation method based 



on the solution of variational inequalities has been used in [3] . 
We start with the formal "interface condition" s = — where we substitute the rhs from (|7k 
with dtu under the smoothness assumption on the solution and with the validity of ([7]) up to the 
interface. Then, we get 

. d'^a{x,t,u) + d^b{x,t,u)+c{x,t,u) 

'^'^ = ' (26) 
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where the rhs is understood in the limit sense x Z' s(t). Again, we consider the transformation 
u = w°' for some a > in ([T]) in the setting ( 25 1 . We obtain 

dtw = ao[dlw + (a - 1)^^] + [bo + 2d,ao]d,w + + -{d^ao + d^w (27) 

w a w a 



where the variables x and t in functions oo, foo and cq are omitted. We are looking for such an a 

)), it 

(28) 



that (27l holds up to the interface x — s{t). Since w{s{t),t) — {'w{x,t) for x s{t)), it 
must hold that 

{a ~ 1)00(9.^)' + ^ for x ^ s{t). 



a 



Due to ao{s{t),t) > S > 0, wc conclude 
dxW — 



Co 



a{a — l)ao 

The negative sign is chosen due to the fact that w is decreasing as x s{t). Since 



m 



dtu dtw 
dxU dxW 



(29) 



(30) 



and \s{t)\ < oo, we choose a so that d^w ^ 0. Then, with respect to (29 1 we choose 

2 

(m-l)a + 2 = =^ oi=- > 1. 

1 — m 



This guarantees that Idxw] < oo. To make use of (|27| in (|30| we have to compute the limit 

2 ^ / (Q-l)ao(9a;w)^ + ^co \ 

xys(t) \ w ) ' 

Both, the nominator and denominator tend to zero as x /" s(i), so the L'Hospital rule can be 
used. We obtain 



(a — l)dxaQ{s{t),t)dx'w + 2{a — l)ao{s{t),t)d'^w + 



1 dxCoisit),t) 
a dxW 



(31) 



{a ~ l)dxao{s{t),t)^ 



2(a- l)aa{s{t),t)dlw 
1 



co{s{t),t) 

a{a — l)ao{s(t),t) a 



dxCo{s{t),t)A - 



a{a — l)aa{s{t),t) 



co{s{t),t) 



Finally, from (26)-(31l we conclude 



3{t) = {2a-l)ao{s{t),t)^ 



a{a- l)ao{s{t),t) ^ 



CQ{s{t),t) 



(32) 



2dxaa{s{t),t) + {a — l)dxCQ{s{t),t) , where a = z— — and u — w" 



co(s(t),i) ■ 



1 — m 



Theorem 3. Consider ([7| under the assumption (25 1. Let the solution u be C^- smooth up to 
the interface s{t). Then, the solution admits an interface s{t) and this is governed by the ODE 
([32]). 

In the case m — 0, the solution u is asymptotically a quadratic polynomial at the interface, 
since u — •uP' , dxW ^ 0, and u is finite on the interface. This is the case for a classical oxygen- 
consumption problem (24 1. 

This result can be extended to the more general case (|7|, with a = aolx, t)g{u)u, b = b^u^ , where 
3(0) > (5 > 0, 7 > 1 and c is as in (pS]). 
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The numerical approximation of the generahzed oxygen-consumption problem (|7| under the con- 
ditions (25 1 proceeds along the same lines as in Section 4, using our governing equation (32 1. The 



approximation d^w on the interface by means of dy 
order polynomial through the points 



is generaly not satisfactory. The third 



(2/Af-3, C'jv-a), (2/JV-2, C'jv-2), (2/JV-1, Cat-i), (1, 0). 



has to be used. 



7 Numerical experiments 

To make an advantage of our numerical approximation a suitable discretization (in the y G (0, 1)- 
variable) can be chosen. Since y — I corresponds to the interface, where the front of the travelling 
and damping semiwave is situated, we increase the density of grid points in the neighbourhood 
of 2/ = 1. At the same time we decrease the density of grid points at the trailing part of the 
damping semiwave type solution to keep the set of all grid points small. This of course depends 
on the solution profile and on the type of the solved model. The distribution of the grid points 
plays an important role in the accuracy and efficiency of the numerical realization. Numerical 
approximation without the interface information requires much more grid points in the space 
discretization then our method with interface modelling to obtain the same order of accuracy. For 
large t the profile of the solution can be significantly changed comparing with that one for small 
t. In that case we can restart the numerical solution of the corresponding ODE-system for large 
t with a new grid (with a changed strategy of grid points distribution) and the corresponding 
initial profile. The main goal is to obtain the same required accuracy with a minimal set of 
grid points, which determines the dimension of our ODE system to be solved. We have used the 
following discretizations. The interval (0, 1) is divided equidistantly in M subintervals and the last 
d intervals are successively subdivided :to 2^, 3^, ...{d — 1)p additional subintervals for p — 1,2, 3. 
We denote these discretizations by Dp. Mostly, we will apply a smooth geometrical dicretization, 
where the length of the (z -I- l)-th subinterval is a g < 1 factor smaller than that of the i-th. At 
this strategy, q is determined from the prescribed number of all grid points N and the length of 
the first subinterval 1/M. We denote this discretization by D4. Space discretization D4 seems to 
be most efficient in most of our numerical applications. 

In our numerical experiments we solve the resulting ODE-system by means of solver odel5s from 
MATLAB hbrary, which is developed for stiff ODE-systems. 



7.1 Barenblatt solution. 

The interface model given by Theorem 1 (the case with cq = 0) coincides with that one in ^. In 
Fig. 1 we present our numerical solution using n = 6 and the discretization with parameters 
M = 10 and N = 20. The corresponding analytical solution cannot be distinguished graphically 
from the numerical one. The time evolution of the interfaces is shown in Fig. 2. The interfaces 
cannot be distinguished, too, even during the long time period (0,200). 

To compare the numerical u„ and analytical Uan solutions, we use a (relative) error estimator 
given by a numerical L2-iiorm . We use discrete time points tj of the solutions approximated in 
Xi = s(tj)yi, i = I, N grid points and define 

{J2f^l\Un{Xi,tj) - Uan{x^,tj)\^{y, - y^^l)s{tj)y/^ 
{E^l \Uan{xi,tJ)\^{y^ ~ y^^l)s{tj)y/^ 

where the space integration is performed over the maximum length of supports of both solutions. 
In the same way we define the Li^rei{tj) error (the exponent 2 is replaced by 1 in the integrands). 
Mostly (in our experiments below) the time evolution of both errors does not differs significantly. 
The time evolution of the relative L2- error is shown in Fig. 3. Using more grid points {N = 



10 



0.7 




X 



Figure 1: Numerical solution of Q with p = 6 in 10 equidistant time sections t S [0, 200] 




time 



Figure 2: Time evolution of the interfaces (numerical and analytical) 
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O 20 40 60 SO 100 120 140 ISO ISO 200 

time 

Figure 3: Time evolution of error 



40, M = 15), the maximum Ls-error is 2.2e - 04 and for [N = 60, M = 20) it is 1.3e - 04. 

Let us approximate Q in a similar way (reducing it to an ODE-system) without the interface 
modelling. We use a mass conservation scheme where 

dlui..,t) ^ ^ ^ ^ [.(.+ 1)C._,.^ + .(^)C|V, - (a(.+ 1) + .(.))Cf )] = C. 

(33) 

with a{i) — Hi — j/i-i, i — I, — 1. At the point j/ = we obtain the ODE 

~ a(l)2 ~ '-^ojL-o 

due to the symmetry arguments {dxu\x=Q = 0). In Fig. 4 we present the corresponding numerical 
solution for p = 6 with an equidistant space discretization of a; € [0, 10] with iV = M = 100 grid 
points. The error L2,rei{t) € (0.13,0.018) for t e (0,200). This error is nearly 100 times higher 
than the error of numerical results based on interface modelling with only = 20 grid points. The 
relative L2.rei(^)"*3rror for the classical numerical approximation increases significantly with the 
reduction of grid points. The most of the contribution to the error is located near the interface. 

In the following Tables 1-3 we present the relation between the space discretization Dp, p — 
1,2,3,4 and the average of L2,rei(^)-error versus the time t E (0,200) in terms of the following 
numerical approximation 

1 

AL — - L2^rel{tj), 

where we take r — 30. We compare our numerical solution (based on interface modelling) with the 
analytical Barenblatt-Pattle solution for n = 6. In Table [T] we present N and AL where M = N 
and discretization Z?4 (i.e., unifom discretization in ?/ g (0,1)). In Table [2] we chose an optimal 
2 < Af < so that AL is minimal. In Table [3] we present N, Dp, d, M and AL where for given N 
and Dp an optimal AI and d are chosen so that AL is minimal. 
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Figure 4: Time evolution of the "classical" numerical solution 
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Table 1: N versus AL Table 2: Optimal N,M versus AL 
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Table 3: Optimal parameters versus AL 



Table 4: N ^ M versus AL 
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7.2 An example of porous media equation with adsorption 



The time evolution of the interface can be very interesting when degenerated diffusion and ad- 
sorption is considered. In the following example we may see the dependence of the time evolution 
of the interface on the local profile of the solution (at the interface) and observe the order of the 
degeneracy very transparently. Let us consider the model problem 

dtu = dluP - Cou'^^P where Cq > 0, 1< p < 2 (34) 

with d^u^lx^Q — 0. For a special initial profile there exists an analytical solution [21] given by the 
formula 

, ' a{tf - x^Yi^ where a{t) = -^-^^——^t +{p- l)a, qr = 

4p-^ p — 1 P ^ 1 

(when the value in parenthesis [.] is nonnegative, otherwise u — 0) which corresponds to the initial 
profile 

[(p-l)a]-9'-(L2-x2)«- ,for|2;|<io 
, for \x\ > Lo 

where a > 0, Lq > are given parameters. 

This example suits to our model considered in Theorem 1 (where n — p,vi — 2 — p, Vp € (1,2)). 
The interface is modelled (see Theorem 1) by 

n 11 1 
s — dxW + Co(n — 1)— — where (3 — , u — , f3 ■ 



uo{x) 



n — 1 OxW P ^ 1 P ~ 1 

Since we know an analytical solution (for the special initial profile above) a very important question 
arises. Does the analytical solution satisfy our interface model given by Theorem 1? 

The answer is positive and this can be verified by some additional computations. This supports our 
results in interface modelling. In Fig. 5-8 we present the analytical solution, numerical solution, 
interfaces and relative error, respectively, for the parameters p — 1.8, N — 100 and M — 80. In 
this case the interface starts to move forwards in the beginning (there is a sufficiently high slope 
of the initial profile at the point Lq — s(0)) and turns backwards at the time t « 5.5, since then 
the influence of the adsorption prevails. 

At time tque G (17, 17.5) the solution is zero and this is a singularity point for our method. 
Also in the neighbourhood oi t = 17 the error significantly increases (the solution is very small). 
In Fig. 9 we present the interfaces using only = 20 grid points. In this case L2.rei{t) € (0, 0.05) 
for t e (0,12.5) and then increases from 0.05 up to 0.15 for t e (12.5,17). As we can see even 
A^ = 20 grid points are sufficient for the numerical approximation. 

Similarly as in the previous example we solve this problem numerically without the interface 
modelling. In Fig. 10 we show the time evolution of the numerical solution when using A^ = 100 
grid points. We consider the domain of the solution to be (0, 20). 

In this case L2,rei{t) G (0,0.05) for t S (0,12) and then increasis from 0.05 up to 0.25 
for t E (12,17). When using only A^ = 40 grid points, the time evolution of the error is 
L2,rei{t) e (0,0.5) for t G (0,8) and further increases from 0.5 up to 5 for t e (8,14). The 
front of the solution is not as sharp as in the case of the Barenblat solution. Here, the best dis- 
cretization strategy is the uniform partition (A^ — M of the interval ye (Oj 1)) • The dependence 
of AL on A^ is presented in Table 4, where the time interval (0, 14) has been considered. 



7.3 Convection-diffusion-adsorption model (Contaminant transport) 

The governing equation of the (simplified ID) contaminant transport model with adsorption is of 
the form 

dtU = Dd^u - dx{vu) - pdtS, in x e (0,L),<> 0, (35) 
14 



Figure 5: The analytical solution in 31 equidistant time sections of t G (0, 16) 
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Figure 7: Time evolution of interfaces (numerical and analytical), N— 100 
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Figure 8: Time evolution of errors, N=100 
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Figure 9: Time evolution of interfaces (numerical and analytical), N=20 




Figure 10: Time evolution of the numerical solution, N=100 
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coupled with the kinetics of adsorption 



dtS = ni'i'iu) ~ S), (36) 

where u is the concentration of the contaminant, v is the velocity of the fluid in porous media 
and S is the adsorbed contaminant by the unit mass of the porous media with the density p. The 
function is an adsorption isotherm, e.g., ^'(s) = cos^,p G (0,1) (Freundlich isotherm). When 



K ^ oo (adsorption is realized in an equilibrium mode), then S — ^'(u) and in this case (351 is 
rewritten in the form 

dtB{u)^ Ddlu~d^{vu), B{u)^u + p^{u) (37) 
which is of the form Q after the transformation (f> = B{u). 

Numerical modelling of the interface (in the case of Freundlich isotherm) has been discussed in 



Section 5. For the construction of the interface model we do not invert B and use (37 1 

B'[u) 

in the interface condition 

dtu 

s{t) = 

In this case the interface is modelled by the ODE 

s{t) = — zdxW for X = s{t), and /? — 



pa(l-p) ' ' l-p' 



u — w^. 



Implementing this interface model and rewritting (37 1 in terms of (j), we proceed as in Section 4 
and obtain the following numerical results. We use the parameters p = 1/2, t G (0,1.9), A'' = 
150, M — 40, V — 1, D — 0.05, a — p—\ and the Dirichlet condition u(0, i) = 1 on the boundary 
X = . The solution (concentration u) is presented in Fig. 11 in 11 equidistant time sections 
(starting from t = 0), where the initial concentration is represented by the first (blue) graph (a 
regularization of the initial zero profile) . The corresponding interface evolution is presented in Fig. 
12. This will be denoted as the case I. Then, the final concentration profile u{x^ 1.9) is taken as an 
initial profile for thecontaminant transport with the homogeneous Dirichlet boundary condition 
m(0, t) = 0, G (0, 40) . The remaining parameters of the model are the same as in the previous 
experiment. This will be denoted as the case II. The evolution of the concentration profile is 
presented in Fig. 13 and the corresponding interface in Fig. 14. 

The numerical difficulty increases when we consider higher order degeneracy represented by 
p = 0.25 and a small diffusion coefficient D = 0.005. We keep other parameters from the previous 
experiment. In Fig. 15-16 we present the time evolution ( in t € (0, 1.9)-case I and t G (0,40)-case 
II) of the concentration. 

The concentration profile drastically changes in Fig. 16. Therefore, we present the corre- 
sponding solution in 11 time sections of the interval t G (0, 10) (in case II) - see Fig. 17. In this 
experiment the considered model is very near to nonlinear transport and the solution is very near 
to the entropy solution of the nonlinear transport with D — Q - see [10] . To increase the density 
of the grid points at the front we take M = 20, N — 150. 

In the fact the result presented in Fig. 17 demonstrates the efficiency of our numerical method. 
The almost piecewise constant initial profile undergoes the nonlinear transport with zero bound- 
ary condition at a: = 0. This shock (at x = 0) is not physically acceptable and immediately 
changes to the rarefaction which moves up to the front shock (physically acceptable) moving with 
the Rankin-Hugoniot speed. The solution endures also after the collision (rarefaction with the 
front shock) . Now we are interested in finding out how many grid points are still sufficient for the 
numerical approximation. In Fig. 18 we present the solution with the same parameters as in Fig. 
17 with only = 60 {M — 10) grid points. The results are nearly the same as the time evolution 
of concentration with N — 200 and M = 20 . This demonstrates the efficiency of our method. The 
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Figure 11: Time evolution of concentration I in 11 equidistant time sections, p:=0.5, N=150 
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Figure 12: Time evolution of interface I, p=0.5, N=150 
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Figure 13: Time evolution of concentrations II, p=0.5, M=60, N=150 
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Figure 14: Time evolution of interface II, p=0.5,M=60, N=150 
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Figure 15: Time evolution of concentration I, p=0.25, D=0.005, N=150 
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Figure 17: Time evolution of concentrations II, p = 0.25, N = 150, D = 0.005, t € (0, 10) 
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Figure 18: Time evolution of concentrations p = 0.25, TV = 60, D = 0.005, t € (0, 10) 
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Figure 19: Time evolution of concentrations a = 1.5, 5 = Ip = 0.95, iV = 100, 1) = 0.005, t G (0, 40) 



experiments with higher degeneracy p < 0.85 are more suitable for our numerical approximations. 
But when p 1 (p = 1 is linear case) troubles arise with our numerical approximations since this 
is the singular case for our method. 

The same mathematical model in ID was treated numerically in |10j using the operator splitting 
method, where in the transport part a semi-analytical solution has been used and in the diffusion 
part the finite volume method has been applied. A very precise numerical solution of this model 
have been applied in a 2D-problem for the transport of contaminant in a dual-well setting. The 
jump in the Dirichlet boundary condition is applied there, to create the pulse shape in the injection 
well and to compute (and also to measure) the response in the extraction well. This gives very 
important information for the solution of the inverse problems (in the determination of the model 
parameters). Since it is an ill-posed problem, a very precise numerical solution is required. There- 
fore our method can be successfully applied there and we will implement it in our forthcoming 
paper. The method used in |10j is limited to Freundlich and Langmuir isotherms. Our method 
also works in a much more general case of sorption isotherms, since it does not depend on their 
geometrical properties. For example, we consider ^'(s) — ^^^^^^ with p E (0,1). This isotherm 
equals to the Freundlich one for 6 = and to the Langmuir one for p = I . In this setting our model 
(37 1 admits the solution with interface {p G (0, 1)). The case p — 1 is reduced to the Langmuir 
sorption isotherm which is a nondegenerated problem. This case is a singular one for our method 
(the speed of interface is oo). However, we take p = 0.95 and expect that our solution will be 
close to the one with Langmuir sorption isotherm solved in |10| with D = (only transport). To 
compare these results, we take D = 0.005, a = 1.5 and b — 1. We use N = 100 and M — 20 in 
Z?4 discretization. The corresponding results are in Fig. 19, which we can compare (graphically) 
with the corresponding results in [10] (Fig. 4). Our concentration profiles at time sections 5 and 
11 then correspond to the ones for < = 16 and i = 40 in [TD]. There is a good agreement. This 
gives us the possibility to obtain a good approximation also in the case of sorption isotherms 
^'(s), 1^*) (0)1 < oo representing a nondegenerated problem (without interface). In that case, in 
the place of ^'(s) we consider an approximation ^'(s^) with p /* 1, which is of course limited by 
the numerical stability of the ODE solver in our setting. We observe that the p — 0.95 leads to a 
good approximation. 

In the next experiment we take b — 0.01 and the remaining parameters are the same as in 
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Figure 20: Time evolution of concentrations a = 1.5,5 = 0.01, p = 0.95, = 100, D = 0.005, t G 
(0,40) 



the previous experiment (with p = 0.95). The resuhs are shown in Fig. 20 and they represent 
a good approximation of hnear sorption (case p = 1). Consequently, when a sorption isotherms 
\E'(s) doesn't leads to the appearance of the interface (nondegenerated case), we can use the 
approximation based on the Langmuir sorption type ^(jip^) with p y 1, b \ close to the 
instability region. 

Now we compare our results with the corresponding ones in |24j (based on the operator splitting 
method) with the following parameters: p = 0.75, D = 0.1, = 1 and a — 2 . We will use the 
discretization parameters TV = 100 and M — 20. Time evolutions of the concentrations are 
presented in Fig 21. The graph corresponding to the 5-th time section can compared with the 
one in [24i (Fig. 1 for t = 15). Then, we use D = 0.001 and other parameters are the same as in 
the previous experiment. Then, we can compare the graph of the solution shown in Fig. 22 (5-th 
time section) with the corresponding results in [53] ( Fig. 2 for t — 15 ). In both cases we obtain 
a very good agreement. 

7.4 Simplified model for turbulent fiow in porous media 

Consider the following mathematical model 

dtU = 92^3/2 _ („3/2 _ ^1/2) 

which is a simplified ID mathematical model for a turbulent flow in porous media -see [16]. Also 
here for a special initial profile Uq, there is an analytical solution given by |16| 

cosh{x/S) 2 ■ , 2(1 y2)e-5*/6 



z(x,t) = y^(tFTT[(l-^F=^=^)+]' with a{t)^ 



^a{ty^ + l (1 + V2)2 - e-5*/3 

where {A)+ = A ioi A > otherwise (A)+ = 0. In this case the interface is given by the formula 

s{t) = 3ln{y/a{t)-^ + l + a-\t)). 

The analytical solution satisfies our interface model given by Theorem 1 (similarly as in Section 
7.2). 
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Figure 21: Time evolution of concentrations a = 2,p = 0.75,7V 100,1? 0.1,7V = 100, M 
20, i e (0,40) 
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Figure 22: Time evolution of concentrations for a — 2,p — 0.75, N — 100,1? — 0.001, iV 
100, M = 20,t e (0,40) 
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Figure 23: Time evolution of the numerical solution for turbulent flow 



Following our approximation method in Sections 2-4, we successively obtain (3 = 2 and the 
ODE-model for the interface evolution (x /' s{t)) 

s{t) = —3dxW — — , where u = up' . 

ZOxW 

We present our results with discretization parameters = 60 and M = 20 in the discretization 
strategy D4. The time evolution of the numerical solution is in Fig. 23. The analytical solution 
graphically coincides with the numerical one in the used scaling. The interfaces are presented in 
Fig. 24 and the errors in Fig. 25. 

As we can see the numerical and analytical interfaces cannot be distinguished in the used 
scaling. The solution is growing since the source term is positive (the term u^^^ dominates u^^"^). 
We also present the solution of a foam drainage model (see [TSl E7] ) 

dtu^dlu^/^ + dx{u'^), n>l. 

We consider the boundary condition dxu\x=o = (symmetry at a; = 0) and the initial condition 
u{x, 0) = uo{x). 

The numerical solution is presented in Fig. 26. 

Finally, we present the numerical solution for the mathematical model of viscous liquid ad- 
vancing over a dry bed. The motion of a thin sheet of viscous liquid over an inclined plate (see 
[5]) is modelled by 

dtu = dlu^ + dx{u^), n>l; 

We consider the boundary condition dxu\x=o = (symmetry at x = 0) and the initial condition 
m(x, 0) = uo(x). 

The numerical solution is presented in Fig. 27. 
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Figure 24: Time evolution of interfaces 




Figure 26: Solution of foam drainage model 
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